Which cells are low quality cells?

Load libraries.

library(Seurat)

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'spatstat.geom':
  method     from
  print.boxx cli 
Attaching SeuratObject
Warning message:
R graphics engine version 14 is not supported by this version of RStudio. The Plots tab will be disabled until a newer version of RStudio is installed. 

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──
✓ ggplot2 3.3.5     ✓ purrr   0.3.4
✓ tibble  3.1.6     ✓ dplyr   1.0.7
✓ tidyr   1.1.4     ✓ stringr 1.4.0
✓ readr   2.1.1     ✓ forcats 0.5.1
── Conflicts ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
x dplyr::filter() masks stats::filter()
x dplyr::lag()    masks stats::lag()

library(pracma)

Attaching package: ‘pracma’

The following object is masked from ‘package:purrr’:

    cross

Load seurat object.

load("/Volumes/Bud_SSD/20210505-seurat/all-aggr-seurat.Robj")

Explore meta data.

head(all.aggr@meta.data)

Plot histogram of number of umi’s.

ggplot(data = all.aggr@meta.data, mapping = aes(x = nCount_RNA)) +
  geom_histogram() +
  geom_vline(xintercept = mean(all.aggr$nCount_RNA), color = "blue") +
  geom_vline(xintercept = median(all.aggr$nCount_RNA), color = "red") +
  scale_x_continuous(trans = "log10")

Looks normal except for bump on left.

Plot histogram of number of umi’s. Color by embryo.

ggplot(data = all.aggr@meta.data, mapping = aes(x = nCount_RNA, color = embryo.id)) +
  geom_histogram() +
  scale_x_continuous(trans = "log10") + 
  facet_wrap(~embryo.id)

Bump on left due to ATU6-45 and BBH6-59.

What are the genotypes of these two embryos? Which experiments did they come from?

print("ATU6-45")

[1] "ATU6-45"

(all.aggr@meta.data %>% filter(embryo.id == "ATU6-45") %>% pull(genotype))[1]

[1] "FIN"

(all.aggr@meta.data %>% filter(embryo.id == "ATU6-45") %>% pull(experiment))[1]

[1] "20171003"

print("AYR4-22")

[1] "AYR4-22"

(all.aggr@meta.data %>% filter(embryo.id == "AYR4-22") %>% pull(genotype))[1]

[1] "FIN"

(all.aggr@meta.data %>% filter(embryo.id == "AYR4-22") %>% pull(experiment))[1]

[1] "20180515"

Both mutant. But from two different experiments.

Plot histogram of number of genes.

ggplot(data = all.aggr@meta.data, mapping = aes(x = nFeature_RNA)) +
    geom_histogram() +
    geom_vline(xintercept = mean(all.aggr$nFeature_RNA), color = "blue") +
  geom_vline(xintercept = median(all.aggr$nFeature_RNA), color = "red")

Looks similar to histogram of number of umi’s. More normal. Except for bump on left again.

Plot histogram of number of genes. Color by embryo.

ggplot(data = all.aggr@meta.data, mapping = aes(x = nFeature_RNA, color = embryo.id)) +
  geom_histogram() +
  facet_wrap(~embryo.id)

Bump is due to same two embryos earlier. ATU6-45 and AYR4-22.

Calculate % mitochondrial genes.

all.aggr[["percent.mt"]] <- PercentageFeatureSet(all.aggr, pattern = "^mt-")
head(all.aggr@meta.data)

Plot histogram of % mitochondrial genes.

ggplot(data = all.aggr@meta.data, mapping = aes(x = percent.mt)) +
  geom_histogram() +
    geom_vline(xintercept = mean(all.aggr$percent.mt), color = "blue") +
  geom_vline(xintercept = median(all.aggr$percent.mt), color = "red") +
  scale_x_continuous(trans = "log10")

Largest portion of distribution normal. Obvious hump of cells with low % mitochondrial genes. Odd that they have low % mitochondrial genes. Usually low is better, but here, most cells have a % between 1 and 10.

Are these cells associated with one embryo?

ggplot(data = all.aggr@meta.data, mapping = aes(x = percent.mt, color = embryo.id)) +
  geom_histogram() +
  scale_x_continuous(trans = "log10") +
  facet_wrap(~embryo.id, scales = "free_y")

Left bump due to more than one embryo. Pervasive feature.

% mitochondrial genes is dependent on number of genes detected. Do these correpond to cells with low numbers of genes?

ggplot(data = all.aggr@meta.data, mapping = aes(x = percent.mt, y = nFeature_RNA)) +
  geom_point(size = 0.1) +
      geom_vline(xintercept = mean(all.aggr$percent.mt), color = "blue") +
  geom_vline(xintercept = median(all.aggr$percent.mt), color = "red") +
  scale_x_continuous(trans = "log10")

So the cells with low mito tend to have lower number of genes. mito seems to distribute normally along number of genes. But if we only filter these out, it doesn’t account for upper range.

Calculate outliers using median absolute deviation.

all.aggr$log.mito <- log10(all.aggr$percent.mt)
median.mito <- median(all.aggr$log.mito)
all.aggr$d.mito <- abs(all.aggr$log.mito - median.mito)
md.mito <- mad(all.aggr$log.mito)
all.aggr$mito.out <- all.aggr$d.mito > md.mito*3
sum(all.aggr$mito.out)

[1] 2256

ggplot(data = all.aggr@meta.data, mapping = aes(x = percent.mt, y = nFeature_RNA, color = mito.out)) +
  geom_point(size = 0.1) +
  scale_x_continuous(trans = "log10")

What does histogram of % mito now look?

ggplot(data = all.aggr@meta.data, mapping = aes(x = percent.mt, color = mito.out)) +
  geom_histogram() +
  scale_x_continuous(trans = "log10")

Genes?

ggplot(data = all.aggr@meta.data, mapping = aes(x = nFeature_RNA, color = mito.out)) +
  geom_histogram()

UMIs?

ggplot(data = all.aggr@meta.data, mapping = aes(x = nCount_RNA, color = mito.out)) +
  geom_histogram() +
  scale_x_continuous(trans = "log10")

Can we find outliers among genes after already filtering out outliers from mito?

median.gene <- median(all.aggr@meta.data %>% filter(mito.out == FALSE) %>% pull(nFeature_RNA))
all.aggr$d.gene <- abs(all.aggr$nFeature_RNA - median.gene)
md.gene <- mad(all.aggr@meta.data %>% filter(mito.out == FALSE) %>% pull(nFeature_RNA))
all.aggr$gene.out <- all.aggr$d.gene > md.gene*3
sum(all.aggr$gene.out)

[1] 221

So 221 cells can be filtered out after already filtering out outliers from mito.

Where do these cells lie along histogram of number of genes?

ggplot(data = all.aggr@meta.data %>% filter(mito.out == FALSE), mapping = aes(x = nFeature_RNA, color = gene.out)) +
  geom_histogram()

Captured those in upper range.

Along number of umis?

ggplot(data = all.aggr@meta.data %>% filter(mito.out == FALSE), mapping = aes(x = nCount_RNA, color = gene.out)) +
  geom_histogram() +
  scale_x_continuous(trans = "log10")

Upper again.

Can we find outliers among umis after already filtering out mito and genes?

all.aggr$log.umi <- log10(all.aggr$nCount_RNA)
median.umi <- median(all.aggr@meta.data %>% filter(mito.out == FALSE) %>% filter(gene.out == FALSE) %>% pull(log.umi))
all.aggr$d.umi <- abs(all.aggr$log.umi - median.umi)
md.umi <- mad(all.aggr@meta.data %>% filter(mito.out == FALSE) %>% filter(gene.out == FALSE) %>% pull(log.umi))
all.aggr$umi.out <- all.aggr$d.umi > md.umi*3
sum(all.aggr$umi.out)

[1] 197

So another 197 cells after already filering out outliers from mito and genes.

Where do these lie along umis?

ggplot(data = all.aggr@meta.data %>% filter(mito.out == FALSE) %>% filter(gene.out == FALSE), mapping = aes(x = nCount_RNA, color = umi.out)) +
  geom_histogram() +
  scale_x_continuous(trans = "log10")

Hard to see where these are.

Plot dots.

ggplot(data = all.aggr@meta.data %>% filter(mito.out == FALSE) %>% filter(gene.out == FALSE), mapping = aes(x = nCount_RNA, y = 1, color = umi.out)) +
  geom_point(position = position_jitter(), size = 0.1) +
  scale_x_continuous(trans = "log10")

So some upper and some lower range cells.

So a total of 197 for umi, 221 for gene, and 2256 for mito totaling 2511 cells. There is some overlap.

But there may be some overlap. Which cells are out for any one of these reasons?

all.aggr$all.out <- (all.aggr$mito.out == TRUE | all.aggr$gene.out == TRUE | all.aggr$umi.out == TRUE)
sum(all.aggr$all.out)

[1] 2511

So those are the cells to filter out.

Where are these cells along mito and genes?

ggplot(data = all.aggr@meta.data, mapping = aes(x = percent.mt, y = nFeature_RNA, color = all.out)) +
  geom_point(size = 0.1) +
  scale_x_continuous(trans = "log10")

ggplot(data = all.aggr@meta.data, mapping = aes(x = percent.mt, y = nFeature_RNA, color = all.out)) +
  geom_point(size = 0.1) +
  scale_x_continuous(trans = "log10") + 
  facet_wrap(~embryo.id)

How about along gene and umi?

ggplot(data = all.aggr@meta.data, mapping = aes(x = nCount_RNA, y = nFeature_RNA, color = all.out)) +
  geom_point(size = 0.1) +
  scale_x_continuous(trans = "log10") +
  facet_wrap(~all.out)

ggplot(data = all.aggr@meta.data, mapping = aes(x = nCount_RNA, y = nFeature_RNA, color = all.out)) +
  geom_point(size = 0.1) +
  scale_x_continuous(trans = "log10") +
  facet_wrap(~embryo.id)

Let’s remove those cells.

all.out <- all.aggr$all.out
save(all.out, file = "/Volumes/Bud_SSD/20210505-seurat/all-out-all-aggr.Robj")

scp /Volumes/Bud_SSD/20210505-seurat/all-out-all-aggr.Robj schea2@hpc3.rcic.uci.edu:/dfs6/pub/schea2/20210505-seurat